Cardiac mechanics Benchmark - Problem 2#
This code implements problem 2 in the Cardiac Mechanic Benchmark paper
Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.
from pathlib import Path
import dolfin
try:
from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
from dolfin import DirichletBC, Constant, interpolate, Mesh
import pulse
import cardiac_geometries as cg
from fenics_plotly import plot
geo_path = Path("geometry")
if not geo_path.is_dir():
cg.create_benchmark_geometry_land15(outdir=geo_path)
geo = cg.geometry.Geometry.from_folder(geo_path)
geometry = pulse.HeartGeometry(
mesh=geo.mesh,
markers=geo.markers,
marker_functions=pulse.MarkerFunctions(vfun=geo.vfun, ffun=geo.ffun),
microstructure=pulse.Microstructure(f0=geo.f0, s0=geo.s0, n0=geo.n0),
)
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 10%] Meshing curve 2 (Line)
Info : [ 10%] Meshing curve 3 (Ellipse)
Info : [ 20%] Meshing curve 4 (Ellipse)
Info : [ 20%] Meshing curve 7 (Ellipse)
Info : [ 30%] Meshing curve 8 (Line)
Info : [ 30%] Meshing curve 9 (Ellipse)
Info : [ 40%] Meshing curve 12 (Circle)
Info : [ 40%] Meshing curve 16 (Circle)
Info : [ 50%] Meshing curve 24 (Ellipse)
Info : [ 50%] Meshing curve 25 (Line)
Info : [ 60%] Meshing curve 26 (Ellipse)
Info : [ 60%] Meshing curve 29 (Circle)
Info : [ 70%] Meshing curve 33 (Circle)
Info : [ 70%] Meshing curve 41 (Ellipse)
Info : [ 80%] Meshing curve 42 (Line)
Info : [ 80%] Meshing curve 43 (Ellipse)
Info : [ 90%] Meshing curve 46 (Circle)
Info : [ 90%] Meshing curve 50 (Circle)
Info : [100%] Meshing curve 63 (Circle)
Info : [100%] Meshing curve 67 (Circle)
Info : Done meshing 1D (Wall 0.0197534s, CPU 0.016493s)
Info : Meshing 2D...
Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info : [ 10%] Meshing surface 13 (Surface, Frontal-Delaunay)
Info : [ 20%] Meshing surface 17 (Surface, Frontal-Delaunay)
Info : [ 20%] Meshing surface 20 (Surface, Frontal-Delaunay)
Info : [ 30%] Meshing surface 21 (Plane, Frontal-Delaunay)
Info : [ 40%] Meshing surface 30 (Surface, Frontal-Delaunay)
Info : [ 40%] Meshing surface 34 (Surface, Frontal-Delaunay)
Info : [ 50%] Meshing surface 37 (Surface, Frontal-Delaunay)
Info : [ 50%] Meshing surface 38 (Plane, Frontal-Delaunay)
Info : [ 60%] Meshing surface 47 (Surface, Frontal-Delaunay)
Info : [ 70%] Meshing surface 51 (Surface, Frontal-Delaunay)
Info : [ 70%] Meshing surface 54 (Surface, Frontal-Delaunay)
Info : [ 80%] Meshing surface 55 (Plane, Frontal-Delaunay)
Info : [ 90%] Meshing surface 64 (Surface, Frontal-Delaunay)
Info : [ 90%] Meshing surface 68 (Surface, Frontal-Delaunay)
Info : [100%] Meshing surface 71 (Surface, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.0214609s, CPU 0.021555s)
Info : Meshing 3D...
Info : 3D Meshing 4 volumes with 1 connected component
Info : Tetrahedrizing 599 nodes...
Info : Done tetrahedrizing 607 nodes (Wall 0.0058539s, CPU 0.005889s)
Info : Reconstructing mesh...
Info : - Creating surface mesh
Info : - Identifying boundary edges
Info : - Recovering boundary
Info : Done reconstructing mesh (Wall 0.0127351s, CPU 0.009543s)
Info : Found volume 2
Info : Found volume 1
Info : Found volume 4
Info : Found volume 3
Info : It. 0 - 0 nodes created - worst tet radius 0.904667 (nodes removed 0 0)
Info : 3D refinement terminated (599 nodes total):
Info : - 0 Delaunay cavities modified for star shapeness
Info : - 0 nodes could not be inserted
Info : - 1934 tetrahedra created in 4.1117e-05 sec. (47036505 tets/s)
Info : Done meshing 3D (Wall 0.0213887s, CPU 0.018004s)
Info : Optimizing mesh...
Info : Optimizing volume 1
Info : Optimization starts (volume = 805.7) with worst = 0.0233132 / average = 0.739813:
Info : 0.00 < quality < 0.10 : 7 elements
Info : 0.10 < quality < 0.20 : 5 elements
Info : 0.20 < quality < 0.30 : 5 elements
Info : 0.30 < quality < 0.40 : 7 elements
Info : 0.40 < quality < 0.50 : 7 elements
Info : 0.50 < quality < 0.60 : 45 elements
Info : 0.60 < quality < 0.70 : 59 elements
Info : 0.70 < quality < 0.80 : 127 elements
Info : 0.80 < quality < 0.90 : 188 elements
Info : 0.90 < quality < 1.00 : 36 elements
Info : 16 edge swaps, 0 node relocations (volume = 805.7): worst = 0.302779 / average = 0.760591 (Wall 0.000340498s, CPU 0.000383s)
Info : No ill-shaped tets in the mesh :-)
Info : 0.00 < quality < 0.10 : 0 elements
Info : 0.10 < quality < 0.20 : 0 elements
Info : 0.20 < quality < 0.30 : 0 elements
Info : 0.30 < quality < 0.40 : 8 elements
Info : 0.40 < quality < 0.50 : 8 elements
Info : 0.50 < quality < 0.60 : 45 elements
Info : 0.60 < quality < 0.70 : 59 elements
Info : 0.70 < quality < 0.80 : 128 elements
Info : 0.80 < quality < 0.90 : 188 elements
Info : 0.90 < quality < 1.00 : 36 elements
Info : Optimizing volume 2
Info : Optimization starts (volume = 805.7) with worst = 0.0233132 / average = 0.736103:
Info : 0.00 < quality < 0.10 : 7 elements
Info : 0.10 < quality < 0.20 : 5 elements
Info : 0.20 < quality < 0.30 : 5 elements
Info : 0.30 < quality < 0.40 : 8 elements
Info : 0.40 < quality < 0.50 : 11 elements
Info : 0.50 < quality < 0.60 : 44 elements
Info : 0.60 < quality < 0.70 : 59 elements
Info : 0.70 < quality < 0.80 : 125 elements
Info : 0.80 < quality < 0.90 : 190 elements
Info : 0.90 < quality < 1.00 : 33 elements
Info : 16 edge swaps, 0 node relocations (volume = 805.7): worst = 0.302779 / average = 0.754569 (Wall 0.000297568s, CPU 0.000308s)
Info : No ill-shaped tets in the mesh :-)
Info : 0.00 < quality < 0.10 : 0 elements
Info : 0.10 < quality < 0.20 : 0 elements
Info : 0.20 < quality < 0.30 : 0 elements
Info : 0.30 < quality < 0.40 : 8 elements
Info : 0.40 < quality < 0.50 : 14 elements
Info : 0.50 < quality < 0.60 : 45 elements
Info : 0.60 < quality < 0.70 : 61 elements
Info : 0.70 < quality < 0.80 : 126 elements
Info : 0.80 < quality < 0.90 : 187 elements
Info : 0.90 < quality < 1.00 : 33 elements
Info : Optimizing volume 3
Info : Optimization starts (volume = 805.689) with worst = 0.0309823 / average = 0.745513:
Info : 0.00 < quality < 0.10 : 7 elements
Info : 0.10 < quality < 0.20 : 3 elements
Info : 0.20 < quality < 0.30 : 2 elements
Info : 0.30 < quality < 0.40 : 6 elements
Info : 0.40 < quality < 0.50 : 11 elements
Info : 0.50 < quality < 0.60 : 42 elements
Info : 0.60 < quality < 0.70 : 61 elements
Info : 0.70 < quality < 0.80 : 128 elements
Info : 0.80 < quality < 0.90 : 188 elements
Info : 0.90 < quality < 1.00 : 33 elements
Info : 12 edge swaps, 0 node relocations (volume = 805.689): worst = 0.319724 / average = 0.76101 (Wall 0.000367529s, CPU 0.000341s)
Info : No ill-shaped tets in the mesh :-)
Info : 0.00 < quality < 0.10 : 0 elements
Info : 0.10 < quality < 0.20 : 0 elements
Info : 0.20 < quality < 0.30 : 0 elements
Info : 0.30 < quality < 0.40 : 6 elements
Info : 0.40 < quality < 0.50 : 11 elements
Info : 0.50 < quality < 0.60 : 42 elements
Info : 0.60 < quality < 0.70 : 61 elements
Info : 0.70 < quality < 0.80 : 132 elements
Info : 0.80 < quality < 0.90 : 187 elements
Info : 0.90 < quality < 1.00 : 32 elements
Info : Optimizing volume 4
Info : Optimization starts (volume = 805.685) with worst = 0.0552945 / average = 0.747184:
Info : 0.00 < quality < 0.10 : 5 elements
Info : 0.10 < quality < 0.20 : 5 elements
Info : 0.20 < quality < 0.30 : 1 elements
Info : 0.30 < quality < 0.40 : 7 elements
Info : 0.40 < quality < 0.50 : 12 elements
Info : 0.50 < quality < 0.60 : 35 elements
Info : 0.60 < quality < 0.70 : 66 elements
Info : 0.70 < quality < 0.80 : 129 elements
Info : 0.80 < quality < 0.90 : 186 elements
Info : 0.90 < quality < 1.00 : 34 elements
Info : 11 edge swaps, 0 node relocations (volume = 805.685): worst = 0.326537 / average = 0.763665 (Wall 0.000487344s, CPU 0s)
Info : No ill-shaped tets in the mesh :-)
Info : 0.00 < quality < 0.10 : 0 elements
Info : 0.10 < quality < 0.20 : 0 elements
Info : 0.20 < quality < 0.30 : 0 elements
Info : 0.30 < quality < 0.40 : 7 elements
Info : 0.40 < quality < 0.50 : 11 elements
Info : 0.50 < quality < 0.60 : 34 elements
Info : 0.60 < quality < 0.70 : 64 elements
Info : 0.70 < quality < 0.80 : 133 elements
Info : 0.80 < quality < 0.90 : 187 elements
Info : 0.90 < quality < 1.00 : 33 elements
Info : Done optimizing mesh (Wall 0.00424632s, CPU 0.002473s)
Info : Optimizing mesh (Netgen)...
Info : Optimizing volume 1
Info : CalcLocalH: 187 Points 472 Elements 370 Surface Elements
Info : Remove Illegal Elements
Info : 169 illegal tets
Info : SplitImprove
Info : badmax = 21.389
Info : 16 splits performed
Info : SwapImprove
Info : 16 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 124 illegal tets
Info : SplitImprove
Info : badmax = 99.5365
Info : 17 splits performed
Info : SwapImprove
Info : 7 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 90 illegal tets
Info : SplitImprove
Info : badmax = 99.5365
Info : 15 splits performed
Info : SwapImprove
Info : 4 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 64 illegal tets
Info : SplitImprove
Info : badmax = 99.5365
Info : 15 splits performed
Info : SwapImprove
Info : 5 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 29 illegal tets
Info : SplitImprove
Info : badmax = 101.081
Info : 8 splits performed
Info : SwapImprove
Info : 0 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 6 illegal tets
Info : SplitImprove
Info : badmax = 101.081
Info : 2 splits performed
Info : SwapImprove
Info : 0 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 0 illegal tets
Info : Volume Optimization
Info : CombineImprove
Info : 19 elements combined
Info : ImproveMesh
Info : Total badness = 1449.12
Info : Total badness = 1358.97
Info : SplitImprove
Info : badmax = 58.0222
Info : 2 splits performed
Info : ImproveMesh
Info : Total badness = 1367.36
Info : Total badness = 1360.57
Info : SwapImprove
Info : 65 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1229.14
Info : Total badness = 1167.69
Info : CombineImprove
Info : 10 elements combined
Info : ImproveMesh
Info : Total badness = 1070.2
Info : Total badness = 1065.73
Info : SplitImprove
Info : badmax = 15.4382
Info : 0 splits performed
Info : ImproveMesh
Info : Total badness = 1065.73
Info : Total badness = 1065.54
Info : SwapImprove
Info : 14 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1051.45
Info : Total badness = 1039.94
Info : CombineImprove
Info : 2 elements combined
Info : ImproveMesh
Info : Total badness = 1021.92
Info : Total badness = 1021.4
Info : SplitImprove
Info : badmax = 15.4382
Info : 0 splits performed
Info : ImproveMesh
Info : Total badness = 1021.4
Info : Total badness = 1021.37
Info : SwapImprove
Info : 7 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1018.11
Info : Total badness = 1013.1
Info : Optimizing volume 2
Info : CalcLocalH: 187 Points 474 Elements 370 Surface Elements
Info : Remove Illegal Elements
Info : 170 illegal tets
Info : SplitImprove
Info : badmax = 21.389
Info : 17 splits performed
Info : SwapImprove
Info : 22 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 119 illegal tets
Info : SplitImprove
Info : badmax = 99.5365
Info : 15 splits performed
Info : SwapImprove
Info : 5 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 85 illegal tets
Info : SplitImprove
Info : badmax = 99.5365
Info : 14 splits performed
Info : SwapImprove
Info : 5 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 54 illegal tets
Info : SplitImprove
Info : badmax = 433.677
Info : 14 splits performed
Info : SwapImprove
Info : 2 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 27 illegal tets
Info : SplitImprove
Info : badmax = 433.677
Info : 9 splits performed
Info : SwapImprove
Info : 1 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 0 illegal tets
Info : Volume Optimization
Info : CombineImprove
Info : 16 elements combined
Info : ImproveMesh
Info : Total badness = 1477.7
Info : Total badness = 1376.18
Info : SplitImprove
Info : badmax = 58.02
Info : 2 splits performed
Info : ImproveMesh
Info : Total badness = 1384.57
Info : Total badness = 1378.78
Info : SwapImprove
Info : 62 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1270.02
Info : Total badness = 1206.98
Info : CombineImprove
Info : 11 elements combined
Info : ImproveMesh
Info : Total badness = 1099.72
Info : Total badness = 1091.5
Info : SplitImprove
Info : badmax = 17.0774
Info : 0 splits performed
Info : ImproveMesh
Info : Total badness = 1091.5
Info : Total badness = 1090.85
Info : SwapImprove
Info : 22 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1068.29
Info : Total badness = 1043.93
Info : CombineImprove
Info : 3 elements combined
Info : ImproveMesh
Info : Total badness = 1020.02
Info : Total badness = 1018.36
Info : SplitImprove
Info : badmax = 15.4382
Info : 0 splits performed
Info : ImproveMesh
Info : Total badness = 1018.36
Info : Total badness = 1018.28
Info : SwapImprove
Info : 7 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1015.73
Info : Total badness = 1012.33
Info : Optimizing volume 3
Info : CalcLocalH: 187 Points 471 Elements 370 Surface Elements
Info : Remove Illegal Elements
Info : 170 illegal tets
Info : SplitImprove
Info : badmax = 14.227
Info : 18 splits performed
Info : SwapImprove
Info : 21 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 121 illegal tets
Info : SplitImprove
Info : badmax = 35.9306
Info : 16 splits performed
Info : SwapImprove
Info : 5 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 88 illegal tets
Info : SplitImprove
Info : badmax = 35.9306
Info : 17 splits performed
Info : SwapImprove
Info : 6 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 49 illegal tets
Info : SplitImprove
Info : badmax = 87.3503
Info : 13 splits performed
Info : SwapImprove
Info : 3 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 16 illegal tets
Info : SplitImprove
Info : badmax = 87.3503
Info : 5 splits performed
Info : SwapImprove
Info : 1 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 3 illegal tets
Info : SplitImprove
Info : badmax = 87.3503
Info : 1 splits performed
Info : SwapImprove
Info : 0 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 0 illegal tets
Info : Volume Optimization
Info : CombineImprove
Info : 23 elements combined
Info : ImproveMesh
Info : Total badness = 1314.97
Info : Total badness = 1256.73
Info : SplitImprove
Info : badmax = 57.7433
Info : 1 splits performed
Info : ImproveMesh
Info : Total badness = 1264.88
Info : Total badness = 1258.71
Info : SwapImprove
Info : 50 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1172.43
Info : Total badness = 1119.49
Info : CombineImprove
Info : 3 elements combined
Info : ImproveMesh
Info : Total badness = 1092.45
Info : Total badness = 1087.72
Info : SplitImprove
Info : badmax = 20.4893
Info : 0 splits performed
Info : ImproveMesh
Info : Total badness = 1087.72
Info : Total badness = 1087.73
Info : SwapImprove
Info : 22 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1053.89
Info : Total badness = 1029.86
Info : CombineImprove
2024-04-11 13:36:06,063 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
Info : 2 elements combined
Info : ImproveMesh
Info : Total badness = 1011.5
Info : Total badness = 1010.13
Info : SplitImprove
Info : badmax = 13.7339
Info : 0 splits performed
Info : ImproveMesh
Info : Total badness = 1010.13
Info : Total badness = 1010.07
Info : SwapImprove
Info : 13 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 997.021
Info : Total badness = 988.982
Info : Optimizing volume 4
Info : CalcLocalH: 187 Points 469 Elements 370 Surface Elements
Info : Remove Illegal Elements
Info : 169 illegal tets
Info : SplitImprove
Info : badmax = 13.4784
Info : 19 splits performed
Info : SwapImprove
Info : 20 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 115 illegal tets
Info : SplitImprove
Info : badmax = 97.9023
Info : 15 splits performed
Info : SwapImprove
Info : 7 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 82 illegal tets
Info : SplitImprove
Info : badmax = 97.9023
Info : 15 splits performed
Info : SwapImprove
Info : 1 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 60 illegal tets
Info : SplitImprove
Info : badmax = 101.081
Info : 14 splits performed
Info : SwapImprove
Info : 1 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 32 illegal tets
Info : SplitImprove
Info : badmax = 193.337
Info : 10 splits performed
Info : SwapImprove
Info : 1 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : 4 illegal tets
Info : SplitImprove
Info : badmax = 193.337
Info : 1 splits performed
Info : SwapImprove
Info : 0 swaps performed
Info : SwapImprove2
Info : 1 swaps performed
Info : 0 illegal tets
Info : Volume Optimization
Info : CombineImprove
Info : 17 elements combined
Info : ImproveMesh
Info : Total badness = 1473.71
Info : Total badness = 1375.72
Info : SplitImprove
Info : badmax = 51.5817
Info : 1 splits performed
Info : ImproveMesh
Info : Total badness = 1379.26
Info : Total badness = 1367.26
Info : SwapImprove
Info : 61 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1274.47
Info : Total badness = 1215.68
Info : CombineImprove
Info : 8 elements combined
Info : ImproveMesh
Info : Total badness = 1133.58
Info : Total badness = 1125.13
Info : SplitImprove
Info : badmax = 21.9556
Info : 1 splits performed
Info : ImproveMesh
Info : Total badness = 1130.48
Info : Total badness = 1129.19
Info : SwapImprove
Info : 21 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 1101.38
Info : Total badness = 1081.27
Info : CombineImprove
Info : 8 elements combined
Info : ImproveMesh
Info : Total badness = 1008
Info : Total badness = 1003.95
Info : SplitImprove
Info : badmax = 12.3103
Info : 0 splits performed
Info : ImproveMesh
Info : Total badness = 1003.95
Info : Total badness = 1003.69
Info : SwapImprove
Info : 13 swaps performed
Info : SwapImprove2
Info : 0 swaps performed
Info : ImproveMesh
Info : Total badness = 996.234
Info : Total badness = 989.321
Info : Done optimizing mesh (Wall 0.19447s, CPU 0.197986s)
Info : 772 nodes 4264 elements
Info : Writing 'geometry/lv_ellipsoid.msh'...
Info : Done writing 'geometry/lv_ellipsoid.msh'
2024-04-11 13:36:06,184 [927] DEBUG UFL_LEGACY: Blocks of each mode:
2024-04-11 13:36:07,024 [927] DEBUG UFL_LEGACY: Blocks of each mode:
2024-04-11 13:36:07,679 [927] WARNING cardiac_geometries.fibers._utils: Unable to write XDMF file:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to find element family.
*** Reason: Element Quadrature not yet supported.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: fd662efc1c0ddefa341c1dac753f717f6b6292a8
*** -------------------------------------------------------------------------
# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["C"] = 10.0
material_parameters["bf"] = 1.0
material_parameters["bfs"] = 1.0
material_parameters["bt"] = 1.0
material = pulse.Guccione(parameters=material_parameters)
# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
return DirichletBC(
V,
Constant((0.0, 0.0, 0.0)),
geometry.ffun,
geometry.markers["BASE"][0],
)
# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
# Solve problem
pulse.iterate.iterate(problem, lvp, 10.0, initial_number_of_steps=200)
2024-04-11 13:36:08,064 [927] INFO pulse.iterate: Iterating to target control (f_166)...
2024-04-11 13:36:08,065 [927] INFO pulse.iterate: Current control: f_166 = 0.000
2024-04-11 13:36:08,065 [927] INFO pulse.iterate: Target: 10.000
2024-04-11 13:36:08,150 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,151 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,151 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,152 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,856 [927] DEBUG UFL_LEGACY: Blocks of each mode:
99 full
2024-04-11 13:36:09,063 [927] DEBUG UFL_LEGACY: Blocks of each mode:
18 full
2024-04-11 13:37:06,773 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,774 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,774 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,775 [927] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,981 [927] DEBUG UFL_LEGACY: Blocks of each mode:
10 full
2024-04-11 13:37:07,170 [927] DEBUG UFL_LEGACY: Blocks of each mode:
3 full
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 182),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 212),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 237),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 262),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 287),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 312),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 337),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 370),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 411),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 460),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 509),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 550),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 583)],
[Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 180),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 210),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 235),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 260),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 285),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 310),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 335),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 368),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 409),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 458),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 507),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 548),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 581)])
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
endo_apex_marker = geometry.markers["ENDOPT"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)
print(
("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
endo_apex_pos[0],
),
)
Get longitudinal position of endocardial apex: -26.5405 mm
epi_apex_marker = geometry.markers["EPIPT"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)
print(
("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
epi_apex_pos[0],
),
)
Get longitudinal position of epicardial apex: -28.1947 mm
V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()